---
title: "7-4 Days alive and free of hospital by day 28"
description: |
This outcome was defined as the number of days alive and free of hospital
from randomisation to 28 days, calculated as 28 minus the number of days
spent in hospital. All patients dying within 28 days will be assigned zero
free days.
author:
- name: James Totterdell
affiliation: University of Sydney
- name: Rob Mahar
affiliation: University of Melbourne
date: today
freeze: true
---
# Preamble
```{r}
#| label: pkgs
#| code-summary: Load packages
library (tidyverse)
library (labelled)
library (kableExtra)
library (cmdstanr)
library (posterior)
library (bayestestR)
library (bayesplot)
library (matrixStats)
library (plotly)
library (lubridate)
library (ggdist)
library (patchwork)
theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ()))
bayesplot_theme_set (theme_set (theme_classic (base_size = 10 , base_family = "Palatino" ) +
theme (panel.grid = element_blank (),
strip.background = element_blank ())))
color_scheme_set ("red" )
options (digits = 4 )
```
```{r}
#| label: data
#| code-summary: Prepare dataset
devtools:: load_all ()
all_dat <- read_all_no_daily ()
acs_itt_dat <- all_dat %>%
filter_acs_itt () %>%
transmute_model_cols_grp_aus_nz ()
acs_itt_nona_dat <- acs_itt_dat %>%
filter (! is.na (out_dafh))
# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
filter_concurrent_intermediate ()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>%
filter (! is.na (out_dafh))
# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
filter_concurrent_std_aspirin ()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>%
filter (! is.na (out_dafh))
# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
filter_concurrent_therapeutic ()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>%
filter (! is.na (out_dafh))
```
```{r}
#| label: models
#| code-summary: Load models
ordmod0 <- cmdstan_model ("stan/ordinal/logistic_cumulative.stan" ) # No epoch or site
ordmod_epoch <- cmdstan_model ("stan/ordinal/logistic_cumulative_epoch.stan" ) # Epoch only
ordmod <- cmdstan_model ("stan/ordinal/logistic_cumulative_epoch_site.stan" ) # Full model
```
```{r}
#| label: helper-functions
#| code-summary: Functions
make_summary_table <- function (dat, format = "html" ) {
tdat <- dat %>%
group_by (CAssignment = factor (CAssignment,
levels = c ("C1" , "C2" , "C3" , "C4" ),
labels = intervention_labels2 ()$ CAssignment[- 1 ])) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafh)),
Deaths = sprintf (
"%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` DAFH, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafh, na.rm = T),
quantile (out_dafh, 0.25 , na.rm = TRUE ),
quantile (out_dafh, 0.75 , na.rm = TRUE ))
) %>%
bind_rows (
dat %>%
group_by (CAssignment = "Overall" ) %>%
summarise (
Patients = n (),
Known = sum (! is.na (out_dafh)),
Deaths = sprintf (
"%i (%.0f%%)" , sum (D28_death, na.rm = TRUE ), 100 * mean (D28_death, na.rm = TRUE )),
` DAFH, Median (Q1, Q3) ` = sprintf (
"%.0f (%.0f, %.0f)" ,
median (out_dafh, na.rm = T),
quantile (out_dafh, 0.25 , na.rm = TRUE ),
quantile (out_dafh, 0.75 , na.rm = TRUE ))
)
) %>%
rename (` Anticoagulation \n intervention ` = CAssignment)
kable (
tdat,
format = format,
align = "lrrrr" ,
booktabs = TRUE ,
linesep = ""
) %>%
kable_styling (
font_size = 9 ,
latex_options = "HOLD_position"
) %>%
row_spec (nrow (tdat), bold = T)
}
make_primary_model_data <- function (
dat,
vars = NULL ,
beta_sd_var = NULL ,
ctr = contr.orthonorm,
p_mult = 2 ) {
X <- make_X_design (dat, vars, ctr)
attX <- attr (X, "contrasts" )
X <- X[, - 1 ]
attr (X, "contrasts" ) <- attX
nXtrt <- sum (grepl ("rand" , colnames (X)))
beta_sd <- c (rep (1 , nXtrt), beta_sd_var)
epoch <- dat$ epoch
M_epoch <- max (dat$ epoch)
region <- dat$ ctry_num
M_region <- max (region)
site <- dat$ site_num
M_site <- max (site)
region_by_site <- dat %>%
dplyr:: count (ctry_num, site_num) %>%
pull (ctry_num)
y_raw <- ordered (dat$ out_dafh)
y_mod <- as.integer (y_raw)
N <- dim (X)[1 ]
K <- dim (X)[2 ]
unique_y <- length (unique (y_mod))
list (N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
M_region = M_region,
M_site = M_site, site = site,
M_epoch = M_epoch, epoch = epoch,
region_by_site = region_by_site,
p_par = p_mult * rep (1 / unique_y, unique_y),
beta_sd = beta_sd)
}
fit_primary_model <- function (
dat,
ctr = contr.orthonorm,
save = FALSE
) {
mdat <- make_primary_model_data (
dat, c ("inelgc3" , "agegte60" , "ctry" ), c (10 , 2.5 , 1 , 1 ), ctr)
snk <- capture.output (
mfit <- ordmod$ sample (
data = mdat,
seed = 813508 ,
chains = 8 ,
parallel_chains = min (8 , parallel:: detectCores ()),
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0 ,
adapt_delta = 0.98
))
if (save) {
mfit$ save_object (file.path ("outputs" , "models" , "secondary" , "7-4.rds" ))
}
mdrws <- as_draws_rvars (mfit$ draws ())
# Add names
epoch_map <- dat %>% dplyr:: count (epoch, epoch_lab)
site_map <- dat %>% dplyr:: count (site_num, site)
names (mdrws$ beta) <- colnames (mdat$ X)
names (mdrws$ gamma_epoch) <- epoch_map$ epoch_lab
names (mdrws$ gamma_site) <- site_map$ site
# Convert to treatment log odds ratios
mdrws$ Acon <- attr (mdat$ X, "contrasts" )$ randA %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))]
mdrws$ trtA <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]
mdrws$ trtC <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c (
"Intermediate vs low" = exp (mdrws$ trtC[1 ]),
"Low with aspirin vs low" = exp (mdrws$ trtC[2 ]),
"Therapeutic vs low" = exp (mdrws$ trtC[3 ]),
"Intermediate vs low with aspirin" = exp (mdrws$ trtC[1 ] - mdrws$ trtC[2 ]),
"Intermediate vs therapeutic" = exp (mdrws$ trtC[1 ] - mdrws$ trtC[3 ]),
"Low with aspirin vs therapeutic" = exp (mdrws$ trtC[2 ] - mdrws$ trtC[3 ])
)
mdrws$ OR <- c (
setNames (exp (mdrws$ trtC), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),
"Ineligible aspirin" = exp (mdrws$ beta[which (grepl ("inelg" , colnames (mdat$ X)))]),
"Age \u2265 60" = exp (mdrws$ beta[which (grepl ("age" , colnames (mdat$ X)))]),
setNames (exp (mdrws$ beta[which (grepl ("ctry" , colnames (mdat$ X)))]), c ("AU/NZ" , "NP" ))
)
return (list (dat = mdat, fit = mfit, drws = mdrws))
}
```
# Outcome Derivation
The outcome is calculated for a patient as:
- missing, if day 28 mortality is unknown (`D28_PatientStatus` ) or if the patient was known to have been alive at day 28 but number of days in hospital (`D28_OutcomeTotalDaysHospitalised` ) is unknown
- 0 if they died by day 28 (`D28_PatientStatus` )
- 28 - min(28, `D28_OutcomeTotalDaysHospitalised` ) otherwise
# Descriptive Analyses
## Overall
As a cross check, @fig-dis-day-vs-days-hospitalised plots the number of days hospitalised (as reported in `D28_OutcomeTotalDaysHospitalised` ) against the day of discharge from the index admission.
```{r}
#| label: fig-dis-day-vs-days-hospitalised
#| code-summary: Compare D28_OutcomeTotalDaysHospitalised to days hospitalised per dischage date.
#| fig-cap: |
#| Number of days hospitalised against day of discharge from index admission.
#| fig-height: 4
pdat <- acs_itt_dat %>%
filter (D28_death != 1 ) %>%
dplyr:: count (D28_day = pmin (28 , D28_OutcomeTotalDaysHospitalised), DIS_day = pmin (28 , DIS_day))
ggplot (pdat, aes (D28_day, DIS_day)) +
geom_point (aes (size = n), shape = 21 ) +
geom_abline () +
scale_x_continuous ("Number of days hospitalised (D28_OutcomeTotalDaysHospitalised)" ,
breaks = seq (0 , 30 , 2 )) +
scale_y_continuous ("Day of discharge from index admission" ,
breaks = seq (0 , 30 , 2 )) +
labs (size = "Count" )
```
The overall distribution of days alive and free of hospital (DAFH) are reported in @fig-dafh-dist for all participants in the ACS-ITT set.
```{r}
#| label: fig-dafh-dist
#| code-summary: Overall distribution of DAFH
#| fig-cap: |
#| Distribution of days alive and free of hospital.
#| fig-subcap: ACS-ITT
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (dafh = ordered (out_dafh, levels = 0 : 27 ), .drop = F) %>%
mutate (p = n / sum (n))
p <- ggplot (pdat, aes (dafh, p)) +
geom_col () +
labs (
x = "Days alive and free of hospital" ,
y = "Proportion"
)
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "outcome-7-4-dist-overall.pdf" ), p, height = 2.5 , width = 6 )
p
```
```{r}
#| label: tbl-7-4-summary
#| code-summary: Summary of DAFH outcome by arm
#| tbl-cap: Summary of days alive and free of hospital to day 28 by treatment group.
save_tex_table (
make_summary_table (acs_itt_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-4-anticoagulation-summary" ))
make_summary_table (acs_itt_dat)
```
## Age
```{r}
#| label: fig-dafh-dist-age
#| code-summary: DAFH by age
#| fig-cap: |
#| Distribution of days alive and free of hospital by age group
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
agegrp = cut (AgeAtEntry, c (18 , seq (25 , 75 , 5 ), 100 ), include.lowest = T, right = F),
dafh = fct_rev (ordered (out_dafh, levels = 0 : 27 )),
.drop = F) %>%
group_by (agegrp) %>%
mutate (p = n / sum (n))
pdat2 <- pdat %>%
group_by (agegrp) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (agegrp, n)) +
geom_col (colour = "grey40" , fill = "grey40" ) +
geom_vline (xintercept = 60 , linetype = 2 ) +
labs (y = "Number of \n participants" ,
x = "Age at entry" ) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 ))
p2 <- ggplot (pdat, aes (agegrp, p)) +
geom_col (aes (fill = dafh)) +
labs (x = "Age" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
geom_vline (xintercept = 8.5 , linetype = 2 ) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-age.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Country
```{r}
#| label: fig-dafh-dist-country
#| code-summary: DAFH by country
#| fig-cap: |
#| Distribution of days alive and free of hospital by country.
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
country,
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (country) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
pdat2 <- pdat %>%
group_by (country) %>%
summarise (n = sum (n))
p1 <- ggplot (pdat2, aes (country, n)) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Country of enrolment" )
p2 <- ggplot (pdat, aes (country, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-country.pdf" ), p, height = 2.5 , width = 6 )
p
```
## Site
```{r}
#| label: fig-dafh-dist-site
#| code-summary: DAFH by site
#| fig-cap: |
#| Distribution of days alive and free of hospital by study site
#| fig-height: 4
pdat <- all_dat %>%
filter_acs_itt () %>%
filter (! is.na (out_dafh)) %>%
dplyr:: count (
Country = factor (PT_CountryName, levels = c ("India" , "Australia" , "Nepal" , "New Zealand" ),
labels = c ("India" , "Australia" , "Nepal" , "New \n Zealand" )),
Site = fct_infreq (Location),
dafh = ordered (out_dafh, levels = 0 : 27 )) %>%
complete (dafh = ordered (0 : 27 ), nesting (Country, Site), fill = list (n = 0 )) %>%
group_by (Country, Site) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup () %>%
mutate (
Country = droplevels (Country),
Site = droplevels (Site)
)
pdat2 <- pdat %>%
group_by (Country, Site) %>%
summarise (n = sum (n)) %>%
ungroup ()
p1 <- ggplot (pdat2, aes (Site, n)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "" ) +
theme (axis.text.x = element_text (angle = 45 , hjust = 1 ),
panel.border = element_rect (fill = NA ))
p2 <- ggplot (pdat, aes (Site, p)) +
facet_grid ( ~ Country, scales = "free_x" , space = "free_x" ) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 1 )) +
theme (axis.text.x = element_text (hjust = 1 , angle = 45 )) +
theme (legend.key.size = unit (0.25 , "lines" ))
p <- p1 / p2 +
plot_layout (guides = 'collect' )
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-country-site.pdf" ), p, height = 4 , width = 6.25 )
p
```
## Calendar Time
```{r}
#| label: fig-dafh-dist-calendar
#| code-summary: DAFH by calendar date
#| fig-cap: |
#| Distribution of days alive and free of hospital by calendar time.
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
yr = year (RandDate), mth = month (RandDate),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (yr, mth) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p1 <- pdat %>%
group_by (yr, mth) %>%
summarise (n = sum (n)) %>%
ggplot (., aes (mth, n)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col () +
labs (
y = "Number of \n participants" ,
x = "Calendar date (month of year)" ) +
scale_x_continuous (breaks = 1 : 12 )
p2 <- ggplot (pdat, aes (mth, p)) +
facet_grid ( ~ yr, drop = T, scales = "free_x" , space = "free" ) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Calendar date (month of year)" , y = "Cumulative \n proportion" ) +
scale_fill_viridis_d ("DAFH" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.25 , "lines" )) +
scale_x_continuous (breaks = 1 : 12 )
p <- p1 | p2
path <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (path, "7-4-calendar-time.pdf" ), p, height = 2 , width = 6 )
p
```
## Anti-coagulation
```{r}
#| label: fig-dafh-dist-anticoagulation
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention.
#| fig-height: 4
pdat <- acs_itt_nona_dat %>%
dplyr:: count (
CAssignment = factor (CAssignment, labels = str_replace (intervention_labels ()$ CAssignment[- 1 ], "<br>" , " \n " )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive.pdf" ), p, height = 3 , width = 6 )
p
```
```{r}
#| label: fig-dafh-dist-anticoagulation2
#| code-summary: DAFH by intervention
#| fig-cap: |
#| Cumulative distribution of days alive and free of hospital by anti-coagulation intervention.
#| fig-height: 4
p <- ggplot (pdat, aes (dafh, cp)) +
geom_step (aes (colour = CAssignment, group = CAssignment)) +
labs (x = "Days alive and free of hospital" , y = "Cumulative proportion" ) +
scale_colour_viridis_d (
"Days alive and \n free of hospital" , option = "A" , begin = 0 , end = 0.8 ) +
guides (fill = guide_legend (reverse = TRUE ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive2.pdf" ), p, height = 3 , width = 6 )
p
```
## Sample Cumulative Logits
Some evidence of a shift from proportional odds at higher values of DAFH.
```{r}
#| code-summary: Plot sample cumulative logits
#| fig-cap: Inspect sample cumulative logits.
trt_counts <- acs_itt_nona_dat %>%
dplyr:: count (CAssignment, out_dafh) %>%
complete (CAssignment, out_dafh, fill = list (n = 0 )) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n))
trt_logit <- trt_counts %>%
group_by (CAssignment) %>%
mutate (clogit = logit (cumsum (p))) %>%
group_by (out_dafh) %>%
mutate (rel_clogit = clogit - mean (clogit)) %>%
filter (out_dafh != 27 )
ggplot (trt_logit, aes (CAssignment, rel_clogit)) +
facet_wrap ( ~ out_dafh) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
ggplot (trt_logit, aes (out_dafh, rel_clogit)) +
facet_wrap ( ~ CAssignment) +
geom_point () +
labs (y = "Relative (to mean) sample cumulative logit" )
```
# Modelling
The primary analysis uses the ACS-ITT set, so all participants have been randomised to the anticoagulation domain and some participants may have been randomised to the antiviral domain.
## Primary Model
```{r}
#| label: fit-model-primary
#| code-summary: Fit primary model
res <- fit_primary_model (acs_itt_nona_dat, save = TRUE )
```
```{r}
#| label: odds-ratio-summary-table-primary-model-save
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters (fixed-effects).
save_tex_table (
odds_ratio_summary_table_rev (res$ drws$ OR, "latex" ),
"outcomes/secondary/7-4-primary-model-acs-itt-summary-table" )
odds_ratio_summary_table_rev (res$ drws$ OR)
```
```{r}
#| code-summary: Odds ratio summary for epoch and site
#| fig-cap: Summary of epoch and site posterior odds ratios.
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" , "7-4-primary-model-epoch-site-terms.pdf" )
ggsave (pth, p, width = 6 , height = 4.5 )
p
```
```{r}
#| code-summary: Posterior contrasts
#| fig-cap: Summaries on comparisons.
plot_or_densities (res$ drws$ compare)
```
:::{.callout-caution collapse="true"}
##### Diagnostics
```{r}
res$ fit$ summary (
c ("alpha" , "beta" , "gamma_epoch" , "gamma_site" , "tau_epoch" , "tau_site" )) %>%
print (n = Inf )
res$ fit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
##### Trace plots
```{r}
mcmc_trace (res$ drws["beta" ])
mcmc_trace (res$ drws["alpha" ])
mcmc_trace (res$ drws["gamma_site" ])
mcmc_trace (res$ drws["gamma_epoch" ])
```
:::
### Posterior Predictive
```{r}
y_raw <- as.integer (levels (res$ dat$ y_raw))
y_ppc <- res$ drws$ y_ppc
y_ppc_raw <- rfun (\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols (acs_itt_nona_dat, tibble (y_ppc = y_ppc_raw)) %>%
mutate (CAssignment = factor (
CAssignment,
levels = names (intervention_labels_short_break ()$ CAssignment),
labels = intervention_labels_short_break ()$ CAssignment))
grp_ppc <- function (grp, d = 0 : 27 ) {
ppc_dat %>%
group_by (grp = {{grp}}) %>%
summarise (
y_lteq = map_dbl (d, ~ mean (out_dafh <= .x)),
ypp_lteq = map (d, ~ rvar_mean (y_ppc <= .x)),
y_eq = map_dbl (d, ~ mean (out_dafh == .x)),
ypp_eq = map (d, ~ rvar_mean (y_ppc == .x))
) %>%
unnest (c (ypp_lteq, ypp_eq)) %>%
mutate (days = d) %>%
ungroup () %>%
pivot_longer (
y_lteq: ypp_eq,
names_to = c ("response" , "event" ),
names_sep = "_" ,
values_to = "posterior" )
}
plot_grp_ppc <- function (dat, lab = "" , xlab = "Probability" ) {
ggplot (dat %>%
filter (response == "ypp" , event == "eq" ),
aes (y = days)) +
facet_wrap ( ~ grp, nrow = 1 , scales = "free_x" ) +
stat_slabinterval (aes (xdist = posterior), fatten_point = 1 ) +
geom_point (data = dat %>% filter (response == "y" , event == "eq" ),
aes (y = days, x = mean (posterior)),
colour = "red" ,
shape = 23 ) +
scale_x_continuous (
xlab, breaks = c (0 , 0.3 , 0.6 ),
sec.axis = sec_axis (
~ . , name = lab, breaks = NULL , labels = NULL )) +
scale_y_continuous ("Days" ) +
theme (strip.text = element_text (size = rel (0.7 )),
axis.title.x = element_text (size = rel (0.7 )),
axis.text.x = element_text (size = rel (0.65 )),
axis.title.y = element_text (size = rel (0.75 )),
axis.title.x.bottom = element_blank ())
}
pp_C <- grp_ppc (CAssignment)
pp_ctry <- grp_ppc (country)
pp_epoch <- grp_ppc (epoch)
pp_site <- grp_ppc (site) %>%
left_join (ppc_dat %>% dplyr:: count (site, country),
by = c ("grp" = "site" ))
p1 <- plot_grp_ppc (pp_C, "Anticoagulation" , "" )
p2 <- plot_grp_ppc (pp_ctry, "Country" , "" )
p3 <- plot_grp_ppc (pp_epoch, "Epoch" , "" )
p4 <- plot_grp_ppc (
pp_site %>% filter (country == "IN" ), "Sites India" , "" )
p5 <- plot_grp_ppc (
pp_site %>% filter (country == "AU" ), "Sites Australia" , "" )
p6 <- plot_grp_ppc (
pp_site %>% filter (country == "NP" ), "Sites Nepal" , "" )
p7 <- plot_grp_ppc (
pp_site %>% filter (country == "NP" ), "Sites New Zealand" , "" )
g1 <- (p1 | p2) / p3
g2 <- p4 / p5 / (p6 | p7)
pth1 <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" ,
"7-4-primary-model-acs-itt-ppc1.pdf" )
pth2 <- file.path (
"outputs" , "figures" , "outcomes" , "secondary" ,
"7-4-primary-model-acs-itt-ppc2.pdf" )
ggsave (pth1, g1, width = 6 , height = 5 , device = cairo_pdf)
ggsave (pth2, g2, width = 6 , height = 6 , device = cairo_pdf)
g1
g2
```
## Sensitivity: Reduced Model (No site term)
The model is re-fit without site terms, but retaining the epoch terms.
```{r}
#| code-summary: Fit reduced model without site term
#| eval: false
fit_no_site_model <- function (dat) {
mdat <- make_primary_model_data (dat, c ("inelgc3" , "agegte60" , "ctry" ), c (10 , 2.5 , 1 , 1 ))
snk <- capture.output (
mfit <- ordmod_epoch$ sample (
data = mdat,
seed = 813508 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0 ,
adapt_delta = 0.98
))
mdrws <- as_draws_rvars (mfit$ draws ())
# Add names
epoch_map <- dat %>% dplyr:: count (epoch, epoch_lab)
names (mdrws$ beta) <- colnames (mdat$ X)
names (mdrws$ gamma_epoch) <- epoch_map$ epoch_lab
# Convert to treatment log odds ratios
mdrws$ Acon <- attr (mdat$ X, "contrasts" )$ randA %**% mdrws$ beta[grepl ("randA[0-9]" , names (mdrws$ beta))]
mdrws$ Ccon <- attr (mdat$ X, "contrasts" )$ randC %**% mdrws$ beta[grepl ("randC[0-9]" , names (mdrws$ beta))]
mdrws$ trtA <- mdrws$ Acon[- 1 ] - mdrws$ Acon[1 ]
mdrws$ trtC <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ OR <- c (
setNames (exp (mdrws$ trtC), c ("Intermediate" , "Low with aspirin" , "Therapeutic" )),
"Ineligible aspirin" = exp (mdrws$ beta[which (grepl ("inelg" , colnames (mdat$ X)))]),
"Age \u2265 60" = exp (mdrws$ beta[which (grepl ("age" , colnames (mdat$ X)))]),
setNames (exp (mdrws$ beta[which (grepl ("ctry" , colnames (mdat$ X)))]), c ("AU/NZ" , "NP" ))
)
return (list (dat = mdat, fit = mfit, drws = mdrws))
}
res <- fit_no_site_model (acs_itt_nona_dat)
```
```{r}
#| code-summary: Odds ratio summary
#| tbl-cap: Posterior summaries for model parameters (fixed-effects).
#| eval: false
odds_ratio_summary_table_rev (res$ drws$ OR)
```
```{r}
#| code-summary: Epoch and site terms
#| fig-cap: Posterior summaries for epoch terms.
#| eval: false
p <- plot_epoch_terms (exp (res$ drws$ gamma_epoch))
p
```
## Sensitivity: Treatment coding
The model is re-fit with priors specified on the treatment contrasts rather than the orthonormal contrasts.
```{r}
#| code-summary: Fit primary model using treatment contrasts
res <- fit_primary_model (acs_itt_nona_dat, ctr = contr.treatment, save = FALSE )
```
```{r}
#| code-summary: Odds ratio table
#| tbl-cap: Posterior summaries for model parameters.
odds_ratio_summary_table_rev (res$ drws$ OR)
```
```{r}
#| code-summary: Epoch and site terms
#| label: odds-ratio-summary-primary-model-epoch-site-terms-prior-sens-trt-ctr
#| fig-cap: |
#| Posterior summaries for epoch and site terms under
#| treatment coding for interventions.
p <- plot_epoch_site_terms (
exp (res$ drws$ gamma_epoch),
exp (res$ drws$ gamma_site),
factor (res$ dat$ region_by_site,
labels = c ("India" , "Australia \n New Zealand" , "Nepal" ))
)
p
```
## Sensitivity: Concurrent controls
Three separate analyses are conducted based on concurrent randomisations as was done for the primary outcome. For these models, the epoch term has been removed.
### Intermediate dose
In the following analyses, the set is restricted to participants who were randomised to either low-dose or intermediate-dose.
- **Set**: ACS-ITT-intermediate
- **Covariates**: anticoagulation intervention, age group, region of enrolment.
```{r}
#| code-summary: Descriptive summary table
save_tex_table (
make_summary_table (acs_itt_concurc2_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-4-anticoagulation-concurrent-intermediate-summary" ))
make_summary_table (acs_itt_concurc2_dat)
```
```{r}
#| code-summary: Fit model
fit_acs_itt_intermediate <- function (fit = TRUE ) {
ctr <- contr.orthonorm
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc2_nona_dat,
contrasts.arg = list (CAssignment = ctr))[, - 1 ]
y_raw <- ordered (acs_itt_concurc2_nona_dat$ out_dafh)
y_mod <- as.integer (y_raw)
N <- dim (X)[1 ]
K <- dim (X)[2 ]
unique_y <- length (unique (y_mod))
mdat <- list (
N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
p_par = 2 * rep (1 / unique_y, unique_y),
beta_sd = c (1 , 2.5 , 1 , 1 )
)
if (fit) {
snk <- capture.output (
mfit <- ordmod0$ sample (
data = mdat,
seed = 75136 ,
chains = 8 ,
parallel_chains = min (8 , parallel:: detectCores ()),
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0
))
mfit$ save_object (file.path ("outputs" , "models" , "secondary" , "7-4-dafh" , "acs-itt-intermediate.rds" ))
} else {
mft <- readRDS (file.path ("outputs" , "models" , "secondary" , "7-4-dafh" , "acs-itt-intermediate.rds" ))
}
mdrws <- as_draws_rvars (mfit$ draws ())
names (mdrws$ beta) <- colnames (X)
mdrws$ Ccon <- as.vector (ctr (2 ) %**% mdrws$ beta[1 ])
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]))
mdrws$ OR <- c (
"Intermediate" = exp (mdrws$ Ctrt),
setNames (exp (mdrws$ beta[2 : 4 ]), c ("Age \u2265 60" , "Australia/New Zealand" , "Nepal" )))
return (list (dat = mdat, fit = mfit, drws = mdrws))
}
res <- fit_acs_itt_intermediate (TRUE )
mdat <- res$ dat
mfit <- res$ fit
mdrws <- res$ drws
```
```{r}
#| code-summary: Sample distribution
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention in subset.
#| fig-height: 4
pdat <- acs_itt_concurc2_nona_dat %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" ),
labels = c ("Low" , "Intermediate" )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh)), colour = NA ) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) +
theme (legend.key.size = unit (0.75 , "lines" ))
p
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-concurrent-intermediate.pdf" ), p, height = 2.5 , width = 4 )
```
```{r}
#| code-summary: Posterior contrast
#| tbl-cap: Posterior densities for comparisons of interest.
p <- plot_or_densities (mdrws$ compare) +
# geom_segment(aes(x = 1, xend = Inf, y = 1 - 0.1, yend = 1 - 0.1),
# arrow = arrow(length = unit(0.5, "lines"))) +
geom_text (aes (x = 1 , y = 1 - 0.1 , label = sprintf ("Pr(OR > 1) = %.2f" , Pr (RV > 1 ))), hjust = 0 , nudge_x = 0.01 ,
family = "Palatino" ) +
scale_x_log10 (breaks = c (0.8 , 0.9 , 1 , 1.1 , 1.2 , 1.3 , 1.4 , 1.5 , 1.6 )) +
labs (y = "Comparison" )
p
```
```{r}
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters.
odds_ratio_summary_table_rev (mdrws$ OR)
save_tex_table (
odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),
"outcomes/secondary/7-4-primary-model-acs-itt-concurrent-intermediate-summary-table" )
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
mfit$ summary (c ("beta" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["alpha" ])
```
:::
### Low-dose with aspirin
- **Set**: ACS-ITT-aspirin
- **Covariates**: anticoagulation intervention, age group, region of enrolment.
```{r}
#| code-summary: Descriptive summary table
save_tex_table (
make_summary_table (acs_itt_concurc3_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-4-anticoagulation-concurrent-stdaspirin-summary" ))
make_summary_table (acs_itt_concurc3_dat)
```
```{r}
#| code-summary: Fit model
ctr <- contr.orthonorm
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc3_nona_dat,
contrasts.arg = list (CAssignment = ctr))[, - 1 ]
y_raw <- ordered (acs_itt_concurc3_nona_dat$ out_dafh)
y_mod <- as.integer (y_raw)
N <- dim (X)[1 ]
K <- dim (X)[2 ]
unique_y <- length (unique (y_mod))
mdat <- list (
N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
p_par = 2 * rep (1 / unique_y, unique_y),
beta_sd = c (1 , 1 , 2.5 , 1 )
)
snk <- capture.output (
mfit <- ordmod0$ sample (
data = mdat,
seed = 135356 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0
))
mdrws <- as_draws_rvars (mfit$ draws ())
names (mdrws$ beta) <- colnames (X)
mdrws$ Ccon <- as.vector (ctr (3 ) %**% mdrws$ beta[1 : 2 ])
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),
"Low with aspirin vs low" = exp (mdrws$ Ctrt[2 ]),
"Intermediate vs low with aspirin" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))
mdrws$ OR <- c ("Intermediate" = exp (mdrws$ Ctrt[1 ]),
"Low with aspirin" = exp (mdrws$ Ctrt[2 ]),
setNames (exp (mdrws$ beta[3 : 4 ]), c ("Age \u2265 60" , "Australia/New Zealand" )))
```
```{r}
#| code-summary: Sample distribution
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention in subset.
#| fig-height: 4
pdat <- acs_itt_concurc3_nona_dat %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" , "C3" ),
labels = c ("Low" , "Intermediate" , "Low \n with aspirin" )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) +
theme (legend.key.size = unit (0.75 , "lines" ))
p
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-concurrent-stdaspirin.pdf" ), p, height = 2.5 , width = 6 )
```
```{r}
#| code-summary: Posterior contrast
#| tbl-cap: Decision summaries on comparisons
plot_or_densities (mdrws$ compare)
```
```{r}
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters.
odds_ratio_summary_table_rev (mdrws$ OR)
save_tex_table (
odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),
"outcomes/secondary/7-4-primary-model-acs-itt-concurrent-stdaspirin-summary-table" )
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
mfit$ summary (c ("beta" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["alpha" ])
```
:::
### Therapeutic dose
- **Set**: ACS-ITT-therapeutic
- **Covariates**: anticoagulation intervention, age group, region of enrolment.
```{r}
#| code-summary: Descriptive summary table
save_tex_table (
make_summary_table (acs_itt_concurc4_dat, "latex" ),
file.path ("outcomes" , "secondary" , "7-4-anticoagulation-concurrent-therapeutic-summary" ))
make_summary_table (acs_itt_concurc4_dat)
```
```{r}
#| code-summary: Fit model
ctr <- contr.orthonorm
X <- model.matrix (
~ CAssignment + agegte60 + ctry,
data = acs_itt_concurc4_nona_dat,
contrasts.arg = list (CAssignment = ctr))[, - 1 ]
y_raw <- ordered (acs_itt_concurc4_nona_dat$ out_dafh)
y_mod <- as.integer (y_raw)
N <- dim (X)[1 ]
K <- dim (X)[2 ]
unique_y <- length (unique (y_mod))
mdat <- list (
N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
p_par = 2 * rep (1 / unique_y, unique_y),
beta_sd = c (1 , 1 , 2.5 , 1 , 1 )
)
snk <- capture.output (
mfit <- ordmod0$ sample (
data = mdat,
seed = 49135 ,
chains = 8 ,
parallel_chains = 8 ,
iter_warmup = 1000 ,
iter_sampling = 2500 ,
refresh = 0
))
mdrws <- as_draws_rvars (mfit$ draws ())
names (mdrws$ beta) <- colnames (X)
mdrws$ Ccon <- as.vector (ctr (3 ) %**% mdrws$ beta[1 : 2 ])
mdrws$ Ctrt <- mdrws$ Ccon[- 1 ] - mdrws$ Ccon[1 ]
mdrws$ compare <- c ("Intermediate vs low" = exp (mdrws$ Ctrt[1 ]),
"Therapeutic vs low" = exp (mdrws$ Ctrt[2 ]),
"Intermediate vs therapeutic" = exp (mdrws$ Ctrt[1 ] - mdrws$ Ctrt[2 ]))
mdrws$ OR <- c ("Intermediate" = exp (mdrws$ Ctrt[1 ]),
"Therapeutic" = exp (mdrws$ Ctrt[2 ]),
setNames (exp (mdrws$ beta[3 : 5 ]), c ("Age \u2265 60" , "India" , "Australia/New Zealand" )))
```
```{r}
#| code-summary: Sample distribution
#| fig-cap: |
#| Distribution of days alive and free of hospital by anti-coagulation intervention in subset.
#| fig-height: 4
pdat <- acs_itt_concurc4_nona_dat %>%
dplyr:: count (
CAssignment = factor (
CAssignment,
levels = c ("C1" , "C2" , "C4" ),
labels = c ("low" , "Intermediate" , "Therapeutic" )),
dafh = ordered (out_dafh, levels = 0 : 27 ),
.drop = F) %>%
group_by (CAssignment) %>%
mutate (p = n / sum (n)) %>%
mutate (cp = cumsum (p)) %>%
ungroup ()
p <- ggplot (pdat, aes (CAssignment, p)) +
geom_col (aes (fill = fct_rev (dafh))) +
labs (x = "Anticoagulation intervention" , y = "Cumulative proportion" ) +
scale_fill_viridis_d ("Days alive and \n free of hospital" , option = "A" , direction = - 1 ) +
guides (fill = guide_legend (reverse = TRUE , ncol = 3 )) +
theme (legend.key.size = unit (0.75 , "lines" ))
pth <- file.path ("outputs" , "figures" , "outcomes" , "secondary" )
ggsave (file.path (pth, "outcome-7-4-descriptive-concurrent-therapeutic.pdf" ), p, height = 2.5 , width = 6 )
p
```
```{r}
#| code-summary: Posterior contrasts
#| tbl-cap: Decision summaries on comparisons
plot_or_densities (mdrws$ compare)
```
```{r}
#| code-summary: Odds ratio summary table
#| tbl-cap: Posterior summaries for model parameters.
save_tex_table (
odds_ratio_summary_table_rev (mdrws$ OR, "latex" ),
"outcomes/secondary/7-4-primary-model-acs-itt-concurrent-therapeutic-summary-table" )
odds_ratio_summary_table_rev (mdrws$ OR)
```
:::{.callout-caution collapse="true"}
#### Diagnostics
```{r}
mfit$ summary (c ("beta" ))
mfit$ diagnostic_summary ()
```
:::
:::{.callout-caution collapse="true"}
#### Trace plots
```{r}
mcmc_trace (mdrws["beta" ])
mcmc_trace (mdrws["alpha" ])
```
:::